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Current aerodynamic designs are often quite complex (geometrically). Flexible 
computational tools are needed for the analysis of a wide range of configurations with 
both internal and external flows. In the past, geometrically dissimilar configurations 
required different analysis codes with different grid topologies in each. The duplicity of 
codes can be avoided with the use of a general multiblock formulation which can handle 
any grid topology. Rather than ‘hard wiring’ the grid topology into the program, it is 
instead dictated by input to the program. 

In the present work the compressible Euler equations, written in a body-fitted finite- 
volume formulation, are solved using a pseudo-time-marching approach. Two upwind 
methods (van Leer’s flux-vector-splitting and Roe’s flux-differencing) have been investi- 
gated. Two types of explicit Solvers (a two-step predictor-corrector and a modified mul- 
tistage Runge-Kutta) have been used with multigrid acceleration to enhance covergence. 

A multiblock strategy is used to allow greater geometric flexibility. The solution 
domain is divided into multiple zones (blocks) and the grid for each block is then 
generated. If the blocks are chosen appropriately, the difficulty of generating a boundary- 
fitted grid is significantly reduced. Also, the placement of wall boundary conditions is 
limited only by the placement of the face of a block. The trade-off for this flexibility 
is the overhead required for the communication between the multiple blocks. In the 
4 present multiblock implementation, two simplifying assumptions have been made. First, 
the interfaces between blocks are assumed to have continuity. Second, the boundary 
condition on any face of a block is assumed to be homogeneous across the entire face 
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(i.e., either completely a wall, an inflow/outflow boundary, or an interface with another 
block). On block faces that have either a wall or an inflow/outflow boundary condition, 
‘standard’ boundary conditions are used. On faces that are interfaces, a special interface 
routine presets the values in two ghost cells (normal to the face) equal to the current 
values in the coincident interior cells in the adjacent block. The updates of the interface 
ghost cells are performed before each iteration in a given block. The iteration of each 
block can then proceed without the need for further information from adjacent blocks. 

There are two possible strategies for the implementation of multigrid with a multiblock 
grid structure: (1) multigrid inside of multiblock and (2) multiblock inside of multigrid. 
With the first strategy, a complete multigrid cycle (or cycles) is (are) performed for a given 
block. Then work begins on the next block and so forth until all the blocks are complete. 
This strategy allows the flexibility of different numbers and/or types of multigrid cycles 
for different blocks, which are adjusted to speed convergence of slowly converging blocks 
(assuming only steady-state results are sought). Unfortunately, communication between 
the blocks is reduced. The rate of convergence is reduced. (The interface boundary 
conditions of the adjacent blocks have to remain fixed in time or a special interface 
condition has to be used for the blocks to communicate from different multigrid levels). 

With the second strategy, multiblock inside of multigrid, all the points on the multigrid 
fine grid (grid h) in all the blocks are updated before the multigrid process continues to the 
coarse grid (grid 2h). There, all the points for all the blocks are updated before proceeding 
to the next coarser multigrid grid. This allows communication between the coarse grids in 
the multigrid cycle. This method can also identically reproduce the convergence history 
of a single block solution using an explicit algorithm- a useful debugging tool for the 
multiblock logic. This latter strategy was used in the present work to compute comer 
flow through a duct, flow from a jet exhaust mixing with freestream air, and transonic 
flow over an ONERA M6 wing, without changing the computer program. 
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SUMMARY 

Several upwind numerical methods for solving the compressible in viscid and viscous flow 
equations are discussed. Due to their explicit nature, the schemes are very simple and easy to apply 
to solutions on multi-block structured grids. Their favourable high frequence damping results in 
good rales of convergence when combined with mnlligrid procedures. Hie schemes are optimized 
using a simple stability and damping factor analysis. Results for three-dimensional test cases are 
shown nnd discussed. Attention is payed to the relative efficiencies of these schemes. 


INTRODUCTION 

After several decades of exislance, computational fluid dynamics finally matured enough to allow 
the solving of complex, realistic problems with confidence. Several numerical schemes are in wide- 
spread use today, solving variety or compressible, viscous and inviscid problems. 

lire simple central difference schemes are particularly popular and effective in inviscid flow 
predictions. The necessity of the addition of ambiguous artificial damping makes them difficult to 
apply to viscous cases. The typical representatives of these numerical methods is for example the 
implicit ARC3D code by Pulliam and Sieger (I) or the proliftc Ilunge-Kutta (R-K) based explicit 
schemes as introduced by Jameson, Schmidt and Turkel (2). 

The upwind schemes are more complex, but also typically more robust and reliable, especially for 
viscous compulations. They are based on several variation of the Riemann problem solver. Typical 
examples are the CFL3D code develo|>cd by Thomas et. al. (.1), based on the Roe's approximate 
Riemann problem solution and Implemented as a flux difference scheme, or the Osher scheme as 
used by Chakrnvnrthy and Osher (1). Roth of these schemes are implicit. The upwind schemes ate 
usually reported to be better suitable for compressible, viscous predictions than the central difference 
rndrs. ORIGINAL PAGE IS 
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The above schemes are very effective In converging lo steady stales on single-block grids of modest 
complexity. Most of lire currently used scheme* are Implicit. The implicit formulation allow* the use 
of large time steps, decreasing the number of llerntions needed lo obtain converged steady slate. Due 
lo the simplification* made during the development of these methods (linearization) and the 
frequent use of explicit boundary conditions, the maximum allowable CFL number has been 
reported as low as ~3.0 for complicated three-dimensional flows. The Implicit part of these solvers 
also provides effective high frequence damping In connection with the application of multigrid 
procedures. 

The application of the above numerical method* lo realistic three-dimensional configurations of 
significant geometric complexity I* usually not possible without the use of multi-block structured 
(zonal) grids. Here, the computational grid I* subdivided Into a number of block* of different size 
and resolution, either overlapping or non-overlapping. A computational grid of this type adapts 
much easier lo the geometric shape of the bodies M well M flow features. The transfer or information 
between the blocks Is typically carried out explicllcly by ensuring the conservation of fluxes accross 
the block interfaces. The consequence of this procedure is a significant reduction of the maximum 
allowable CFL number lo values between 2. and 6. At these CFL numbers, an explicit upwind 
scheme of good accuracy when applied to viscous problems and suitable for multigrid procedures 
seems to offer a belter choice. 

'I here is a very large number of explicit schemes that have been used or are still in use for solving 
the compressible flow equations in their Invisrid form (Filler) or viscous form (Navicr-Slokes). The 
numerical method should be simple, efficient, robust, have effective damping of high frequence errors 
(ncccesary for mulligrid), have low dispersion (low phase error will reduce spurious oscillations and 
result in faster convergence fates), low levels of numerical dissipation for accurate predictions of 
viscous effect* and maintain high resolution on stretched grids. 


II ASIC ALGORITHM 

The scheme is based on Floe's flux difference splitting (5,3). This type of upwind scheme is 
simplified by linearizing the tllcmann problem between two cell interfaces about a slate obtained by 
Hoe s averaging procedure, fief. (5). 1 lie present discretization employed the finite volume approach, 
with the slate variables at the cell Interface* determined by the MUSCL interpolation, using mostly 
thefo-called k scheme. 

In the in viscid case, the governing equations are the three-dimensional Euler equations, 
expressed in nondimenslonal conservation law form for curvilinear coordinates (, tj, < , using the 
usual Huid dynamics notation, as; 
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where Q Is Hie veclor of dependent conservative flow variables and F, (J and fl Arc the usual flux 
vectors in the direction (, rj And respect. The semidtscrele form of the governing equations is 
obtained ffom (1) by carrying out the spatial discretization, assuming = 1.0 

^( wF Mj' F l| (2) 

In equation (2), the numerical flux vector at the cell interface between cells lj,k and 

i + lj,k; ij«k tlie spatial indices In the directions (, r; and respectively. The expressions in the 

other two directions arc similar and will not be repented here. Then, 

t= n rw= , ^|- r ') ,+, %i c! ii ,+<, v.r n M> < s > 

Following fief. |5j f the numerical flux | J ’ f° f Is evaluated ns 

F i4) = H F a + F r. S (l A (l s i , ( t »n c Jr.)l W 

Here, = F(Qj^) * = F(Q| ) » *here and Qj nrc the conservative (low variables 

describing the slates right and left of the cell Interface and are obtained by an appropriate 
inerpolation or extrapolation. "Ihe last expression In eq. (4) is a dampimg term due to the upwind 
character of the present approach; and A^ are the eigenvector and eigenvalue matrices of the flux 
Jacobian matrix evaluated using the Hoe averaged flow variables p, 0, * # ft and fi, given in (5). 

The selection of a particular type of time stepping will determine the characteristics of the 
numerical method. In the present work, explicit multi-stage schemes were considered; their S stage 
version can be written In a somewhat general form as 



1 he above formulation of the time stepping procedure applies not only tb the modified Runge- 
Kulta method in (2), but also many other existing schemes. From these infinitely many possible 
schemes, only a few are subside for Ihe aIk>vc flux difference splitting spatial discretization. They 
were selected by considering the much simplified scalar, linear case of the wave equations. 
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In order lo compare Rome of Ibe time-stepping procedure* (B), the standard Fourier stability 
analysis of several scheme* was carried out, yielding the amplification factor g and the phase error 
for the simple wave equation 

$? + = 0 ■ '>• w 

There are two major block* of scheme variations that fit eq. (5). 

Modified Beam-Warming f DVV) Scheme 

This acheme wa.s originally introduced hy Hearn and Warming [fij. It ia a two-step method 
consisting of a first order predictor and second order corrector. It ha.* been slightly modified lo make 
it suitable for the present finite volume formulation: 

Predictor: 

«!, = «" i ( >h - i (»*) 

q . 1 = q." . n* 

Corrector: 

<*l = Q" + l(Qi' • «"|) i Qit = Q"+i - i(«"+ 2 • Q/+i> <">> 

q . 2 = q" +i = qj 1 - ai n 2 

In both steps, only tire extrapolations (7a) and (7b) are different; the algebra that completes one 
iteration, consisting of eq. (3) with (4), is identical. The stability analysis of this scheme is given in 
every good textbook on CFD and will not be repealed here. The resulting plot of tire magnitute of 
the amplification factor |g| can be seen In Fig. la, plotted a* a function of tbc spatial wave number 
0 between 0 and n and the CFb-number <r between 0 and 2. 

The scheme is still second order accurate in apace and lime, satisfies the "shift condition" and has 
a stability limit of <r = 2. It has a relatively modest disspertion. It was named the "1-2 BW scheme". 
It is a very simple ami efficient method, consisting of only two steps. It performed well on single 
grids, in particular when applied to unsteady flow predictions. 

There were, however, some problems with this scheme. Since the two steps are different, the 
steady stale result will depend on the time step. This phenomen was observed tn only a few cases, 
represented by convergence lo residual that was larger than "machine mo". At this time, none of 
the flux limilors tested In this scheme converged more than two orders of magnitude. The resulting 
flow fields agreed well with other, fully converged numerical results. Thin type of behaviour has been 
observed by oilier Investigators, but Is still disturbing. 

A more serious problem Is the Increase of the damping factor lo 1.0 at high frequency and <r = 1. In 
4 scafar case, tbc CFL number can be kept at Its optimum value (1.7 in the case of the 1-2 BW 
scheme), but In the ease of the Euler equations there are three dislinc eigenvalues in each direction. 
Typically, local time stepping will be Implemented, where each cell will he advanced at its optimum 
lime step, corresponding lo the maximum eigenvalue at that cell. This means that only (lie 
maximum eigenvalue will correspond to the optimum CFL number for high frequence damping and 
one or more might rorre*pnnd to the CFF, number ranges willi vary liltlr damping. Tins was 
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manifested by the Inc V of convergence of the 1-2 HW scheme when utilhed in a Multi-Grid 
procedure. For mtilligrid applications, this scheme was modified by making the first step (predictor) 
also second order accurate: 

q|, = Q" + KQf • Of,) i 'lit = Q"n • l(Q" + 2 ' 1 1) (») 

'I lie plot of the damping characteristics of this scheme, cnllrd 2 2 I1VV, Is shown in Fig. lb. The 
maximum stable CFb number 1* now only 1.0, but the high frequence damping i s much belter. This 
scheme worked well with Multi-Grid and will be discussed below. 


Modified Rungc- Kiitt* Methods 

'Ihe modified R-K methods, with the standard set of coefficients, have been rather successful! in 
combination with the central difference spatial discretization. They have been, however, performing 
very poorly with upwind differencing. 'Ihe standard roefncienls hnve to be modified to achieve 
belter performance, resulting In schemes that are In general of reduced nccurary in lime. In order to 
find the optimum sets of the R-K coefficients, a Fourier stability analysis wiw carried out, similar to 
Hie approach in (2). 'Ihe extrapolation of the stale variables to the cell interfaces was based on the 
so called a-sebrme: 


Qf, = Qf ’ + M(l - + (1 b ; (Da) 

Aj = Q,-Q M ; A. + = Q. +I -Q i; (9b) 

Qri = Q i+l • \'IM' + * f i+l> A i+l + 0 - M 

where l is one of the possible flux lirnitors. The value of the parameter k determines the spatial 
accuracy of the scheme; is fully upwind, second order accurate; k= 0 is the upwind biased 

second order Fromm scheme; is upwind biased Hurd order and rc=l is second order central 
difference scheme, lire first order scheme Is obtained by setting £=0. In the present stability 
analysis, the limilor l was set to I (no limifor) for simplicity. 

The It- K time stepping for equation (6) can l>e written as 

u" = n n - n* ‘ j ii" = (io) 

llrrr, n<,=1 for consistanry, wllb ii”* 1 = u^. Making the assumption Ibnl u in harmonic, leads lo 

tis '4x = fl = "o*- 71 = *«* (H) 

Defining P=.sl *, substituting eq. (I I) In eq. (10) and some simple algebra yields 


K = (<l ^ /u n ) = I 4- P + T ••• k OjOj-.-Og |I 


,S 


( 12 ) 

Clearly, for stability, 9t(P)<0 and |g|<l. The stability and damping properties of the scheme are 
associated with the complex polynomial (12). T he damping |g| is a function of Hie complex P=xd ly, 
l=jn, and can be best shown as a plot of contours of constant |g| between 0 and 1. However, 

"I = = «.' W (13) 

where ft is the spatial wave number, ranging from 0 to rr. Values of ft between | and x are 
considered high frequencies. F now represents the Fourier transform of the spatial differencing 
operator and can be superimposed on the damping factor contour plots. Defining 


»=w 
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( 15 ) 


m llie CFL number, and observing tlml | =cu |^| f* vc * 

A = -m = "j.j) 

Eq. (15), combined with eq.(I3) and (9), for f=l, yields finally 

P = -rr(l-e l/? ){l + ( 1 C) 

P Is a function of <r and fl\ the plot of P(<r,/?) and fg| can he used to optimize the coefficients a 8 . 
The resulting stability plots arc shown only in the tipper half of the negative real part of P (fourth 
quadrant) since they are symmetric with respect to the y =0 coordinate line. The optimization of the 
coefficients n n was carried out by displaying these stability plots on a PC-type microcomputer. The 
changes in the shape of the |g| contours were observed in real time as the coefficients were changed. 
The "islands" of the low value of |g| correspond to the roots of the polynomial (12). The main 
purpose of this optimization was to find a combination of the coefficients cv 8 such that, for as large 
<r as possible, there would he good high frequency damping (low value of |g|) for a large range or <r 
(maximum size of the "islands" as close to the real axis as possible). The optimization was 
performed for four different spatial discretizations (l-*t order; 2 -nd order fully upwind, «=d; 2 nd 
ordrr Fromm scheme, *= 0 ; ,1 rd ordrr upwind bias, K~j|) for the two-, three- and four-stage R-K 
schemes. Even in the case of the four stage scheme, the actual optimization of the three coefficients 
(crg=l) was relatively easy and quick, once the influence of the coefficients on the stability plots was 
understood. Only a few selected stability plols of the most interesting schemes will he presented 
here. The optimized coefficients are summarized in Table I. 

1 he simplest schemes to optimize were all the two singe versions, since only one coefficient is 
freely selectable. The case of *=0 is shown In Fig. 2 a. For the optimum coefficient o-j = 0.-12, the 
theoretically determined maximum slahle cfl number was ]. The Fromm scheme performed very 
well in most of the lest cases due to Its very low numerical dispersion. 

The four-stage R-K schemes are more Interesting. Here, the optimum combination of three 
coefficients lias to he found. The standard R-K coefficients (r>| = |; a j — a 3 = |» o 4 = 1 )» shown in 
Fig. 2b for the third-order scheme, performeml, as expected, relatively poorly. H should be noted 
that, when multi grid procedures are Implemented, the maximum CFL-number is of little 
importance. The high frequency damping (or lack of it) will effect the rate of convergence to steady 
stale more significantly than the CFL-nmnber. In this case of standard coefficients, the maximum 
stalde CFL number Is relatively low at <z = l.7, hut, more significantly, the damping of high 
frequencies of error propagation (at 0>%) Is very weak. The result of the present optimization, 
shown in Fig. J l, is much more promissing. 

It should be mentioned here that, in a parallel effort, van Leer et. al. ( 7 ) also tried to optimize the 
R-K coefficients for applications with the upwind methods. Their approach was somewhat different: 
assuming that a genuine and practical multi dimensional characteristic formulation of the Euler 
equations could he found, they optimized the R-K coefficients for only one value of rr, argueing that 
each wave would he prnpagalrd at its optimum CFI^nntnl**r. Unfortunately, there no such 
formulation for the fhrrr dlmrnnonat ra«*\ '! h* advwtifaf** of flu* approa'h waa thnf flir rrlrrlion 
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process could be Automated, eji^dnatlng the ’’guessing game” involved In the present approach, 
1 heir results are shown In Fig. 1W and X Generally, their maximum tt was lower, and the damping 
effective over a narrower range of CFL values. ^ 

In the case of the third order scheme, the present optimum rt-s are shown in Fig.'ll In the real 
test cases, however, the van Alhada lirnilors were implemented, shifting the P-curvc to the left. The 

real optimum CFL number was therefore much lower (1.8). A new set of coefficients for the limited 

*9 

case was found. It is shown in Fig.^ (o|=.ll; o ? =. 215; o 3 =. / I8); the maximum <r was 2.5. 


RESULTS 

Most of the above schemes were tested on number of two* and three-dimensional cases. In general, 
Ibe optimum CFL numbers agreed well with the above theory. In supersonic rases, the real optimum 
CFL number was usually somewhat higher; In transonic eases, it was mostly the same. 

Due to the space limitations of this paper, only one three-dimensional case will be discussed. Here, 
the supersonic flow in a channel with two 10° compression ramps, at the bottom anti left vertical 
walls, with an Inflow Mach number of ,1.17, formed a conical shock flow with two Mach tripple 
points. The twilling flow field, obtained by the Fromm scheme, is shown in Fig. 3 . The 
compulations were carried out using a 32x12x32 body filled grid. The Full-Mulli-Grid procedure 
(FMG) was used with three grid levels with two iterations on each grid level. Considering the results 
of the above stability analysis, the most promissing schemes of practical importance were the four- 
stage Fromm scheme (*=0) and the third order biased scheme (a=J). The Fromm scheme was 
preferred due to its low numerical dispersion, demonstrated by resulfs with the least oscillations (no 
lirnilors). The convergence of the k=0 scheme with r> f =.ll, ojrr.255, o 3 =.4fl for different CFL- 
numhers Is shown In Fig. 4. The best rate of convergence was achieved at CFL numbers between 2.0 
and 2.3, which wm in a surprising agreement with the theory. This scheme, at a (7=2.0, is compared 
with the simple 2-2 RW scheme at Us maximum <7 = 1.0 in Fig. 5. The 2-2 RW scheme seems to be 
much slower, but it is only a two stage nr. Iienie and thus needs only about half the CPU time per 
one iteration. Consequently, it Is approximately 50% less efficient. 

Finally, the present *=0 scheme is compared with the optimised version of the *=$ scheme in 
Ref. (7), used at the optimum CFL number of 1.732. Its performance is much worse, since only one 
of the three eigenvalues corresponds to the optimum damping case. The application of the van 
Alhada limllor lead to a limit cycle with no apparent convergence. This problem has been reported 
by other Investigators but was still disapolnling, since all* the two-dimensional cases converged with 
the limitor In effect. The problem is being further Investigated. 


CONCLUSIONS 

The present study Investigated several types of explicit upwind schemes for solving the 
rnnipre**ih|e flow problem*. The simple 17 RW scheme seems to very elTrrlive for predicting 
unsteady flows, foil U dor* not work with Multi fJtbl due to its Imuffirlenl damping of high 
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frequencies at wide range of CFL numbers. The most promising schemes for Multi-Grid 
computations are multi stage R-K methods with optimised sets of coefficients. Here, versions with 
coefficients optimized for wider range of CFL numbers seem to be more robust than specialized 
versions, lire k ~ 0 and a = schemes were successfully tested on several two-dimensional cacs and 
subsequently included in a multi-block, multi grid code used for predictions or complex three- 
dimensional flows. 
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